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Abstract 



We develop a fourth order simulation algorithm for solving the stochastic 

00 ■ 

Langevin equation. The method consists of identifying solvable operators in 
C the Fokker-Planck equation, factorizing the evolution operator for small time 

o 

steps to fourth order and implementing the factorization process numerically. 

A key contribution of this work is to show how certain double commutators in 

the factorization process can be simulated in practice. The method is general, 

applicable to the multivariable case, and systematic, with known procedures 
• i— i . 

for doing fourth order factorizations. The fourth order convergence of the 



resulting algorithm allowed very large time steps to be used. In simulating the 
Brownian dynamics of 121 Yukawa particles in two dimensions, the converged 
result of a first order algorithm can be obtained by using time steps 50 times as 
large. To further demostrate the versatility of our method, we derive two new 
classes of fourth order algorithms for solving the simpler Kramers equation 
without requiring the derivative of the force. The convergence of many fourth 
order algorithms for solving this equation are compared. 
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I. INTRODUCTION 



A stochastic differential equation of the form 

x i = Gi(x) + D ij Z j (t), (1) 

or its equivalent Fokker-Planck equation 

^P(x,t) = LP(x,t) = [^Dijdidj - a i G 4 (x)]P(x,i), (2) 

is used to describe a variety of physical and chemical processes Jl|. Even in the Langevin 
case, where the diffusion matrix is position independent, it is difficult to derive numerical 
algorithms for solving it beyond second order . A direct Taylor expansion || approach is 
laborious, giving no insight into the overall structure of the algorithm and requires an eight 
term expansion to achieve 4th order accuracy Heretofore, no fourth order Langevin 
algorithm has been derived and applied to systems of more than one particle. 
The Fokker-Planck equation (Eh can be formally integrated to give 



iV 



P(x, t) = e* L P(x, 0) = [e eL J P(x, 0). (3) 

This equation can be solved by factorizing the short time Fokker-Planck evolution operator 
e eL = e e ( T+D ^ into exactly solvable parts. In this work, we will take = Sij and define 
operators 

T= X -dA and D = -diGi(x), (4) 

with implied summations. This idea of operator factorization is not new, and has been 
used to derive a number of second order Langevin algorithms [|5]||. We will briefly review 
the basic idea in Section II. However, it is only recently that one learns how to factorize 
operators of the form e t ( T+D '> to fourth order with positive coefficients |||9| . All such fourth 
order factorizations require the evaluation of the double commutator [D, [T,D]], which is 
rather formidable at first sight. We will show in Section III, how this commutator can be 
implemented judiciously to yield a fourth order Langevin algorithm. To demonstrate the 



high order convergence of this algorithm, we use it to simulate the Brownian dynamics of 121 
Yukawa particles in two dimensions, a system that has been studied extensively by Branka 
and Heyes [|1(J using second order algorithms. 

To further demonstrate the utility of the factorization method for solving stochastic 
equations, we derive systematically a number of fourth order algorithms for solving the 
Kramers equation in Section IV. Drozdov and Brey fll] have used a similar factorization 
method to solve this equation in one dimension using grid points. Hershkovitz Q has also 
derived a fourth order algorithm by Taylor expansion. In both cases, it is not obvious how 
their respective approaches can be generalized to the multivariable case. We give a detail 
comparison of all algorithms using Monte Carlo simulation, which can be easily generalized 
to any dimension. Finally, we summarize our findings and present some conclusions in 
Section V. 



II. OPERATOR FACTORIZATION 

When the operator e eT acts on P(x, t), it evolves the latter forward in time according to 
the diffusion equation 

^P(x,t) = ^P(x,t). (5) 

If {xi} is a set of points distributed according to P(x, t), then the distribution e time later 
can be exactly simulated by updating each point according to 

x[ = x i + v / e&, (6) 

where {^} is a set of Gaussian distributed random numbers with zero mean and unit vari- 
ance. When the operator e eD acts on P(x, t), it evolves the latter forward in time according 
to the continuity equation 

£p(x,f) = -a[G 4 (x)P(x,f)], (7) 
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where Gj(x)P(x, t) = Jj(x) is the probability current density with velocity field Gj(x). The 
continuity equation can also be exactly simulated by setting 

< = (8) 

where Xj(e) is the exact trajectory determined by 

f = G(x), (9) 

with initial condition Xi(0) = Xj. 

Thus, if e e( - T+D * > can be factorized into products of operators e eT and e eD , then each such 
factorization will give rise to an algorithm for evolving the system forward for time e. For 
example, the second order factorization, 

e h T e tD e^ T = exp[e(T + D) + 0(e 3 ) ■ • •], (10) 

leads to a second order Langevin algorithm || 



Vi = Xi + ^e/2, 

x' 4 = 2 / i (e)+C , v / ^2, (11) 

where £j and ^ are independent sets of zero mean, unit variance Gaussian random numbers. 
For a second order algorithm, it is sufficient to solve for the trajectory yi(e) correctly to 
second order in e, e.g. via a second order Runge-Kutta algorithm: 

yi (e)=y l + eG l (y+ l -eG(y)y (12) 

Alternatively, one has the factorization, 

e ^V T e^ D = exp[e(T + D) + 0(e 3 ) • • •], (13) 

which yields the second order algorithm 

Vi = Xi(e/2) +&v / e> 

< = y l (e/2). (14) 



Again, it is sufficient to solve the trajectory equations Xi(e/2) and yi(e/2) correctly to second 
order via the Runge-Kutta algorithm. Despite the appearance that this algorithm requires 
solving the trajectory equation (|^) twice, it can be shown || that by expanding the two 
trajectories to second order and recollecting terms, one arrives at the second order Runge- 



Kutta Langevin algorithm However, the canonical form of (13), with two evaluations 

of the trajectory, usually has a much smaller second order error coefficient. 

The method of operator factorization thus appears to provide a systematical way of 
generating higher order algorithms. Unfortunately, Suzuki |12| proved in 1991 that, beyond 
second order, for any two operators, T and D, it is impossible to factorize the evolution 
operator as 

N 

exp[e(T + £>)] = II exp^eT] exp[6 4 eL>] (15) 
i=i 

for any finite N, without having some coefficients and hi being negative. In the present 
context, since e aieT is the diffusion kernel, a negative a« would imply that one must simulate 
the diffusion process backward in time, which is impossible. Thus factorizations of the form 
(|T5]) cannot be used to derive higher order Langevin algorithms. 



III. A FOURTH ORDER LANGEVIN ALGORITHM 

The essence of Suzuki's proof is to note that in order to obtain a fourth order algorithm, 
one must eliminate third order error terms involving double commutators [T, [Z?,T]] and 
[D, [T,D]}. With purely positive coefficients a, and 6j, one can eliminate either one or the 
other, but not both. Thus to obtain a fourth order factorization with all positive coefficients, 
one must retain one of the two double commutators. Recently, Chin has derived three 
such factorization schemes, two of which were also found previously by Suzuki ||. 

The form of the operators T and D, as given in ([|), dictates that one should keep only 
the commutator [D, [T, D}}, which is at most a second order differential operator. Since the 
velocity (or force) field G is usually given in terms of a potential function V(x), 
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Gi(x) = -diV(x) 



(16) 



the double commutator has the form 



[D, [T,D]] = a i d j f iJ + dM 



(17) 



where 



fi,j = Vi,j,kVk — 2Vi y kVj y k 

1 



(2Vi,j,kVj,k + Vi,jVj,k,k — Vi,j,k,kVj) ■ 



(18) 



The indices on V indicate corresponding partial derivatives. Since the operator D requires 
solving for the particle's trajectory, we must minimize its occurrence. This dictates that we 
use a variant of Chin's scheme B 101 to factorize 



exp [e (T + D)] = exp 



exp | -D) exp Uf 



x exp ( —D ) exp 



+ 0(e 5 ), 



where we have included the double commutator in T 



(19) 



f = T+-(2v / 3-3)[A [T,D]]. 



To obtain a fourth order algorithm, we must simulate this new term 

<r3 



exp i±f) = exp 



-j=T + € - (2 - V3) («Wm + ft 



(20) 



(21) 



V3 24 

correctly to 4th order. If we simply took all x dependent terms in this operator as fixed, 
evaluated at the starting point, this operator would describe a non- uniform Gaussian random 
walk. However, this normal ordering would be correct only to third order. To implement it 
to fourth order, we first decompose it as 



exp (^) =exp fe 



T exp 



6 

24 



(2 - v^) {didjfij + d { 



exp 



2v/3 



T +0(e 5 ). 



(22) 
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If fij is positive definite, normal ordering the middle operator above, i.e. interpreting it 
as a non-uniform Gaussian random walk with fij evaluated at the starting point, would 
be correct to 4th order (actually to 5th order). However, if some eigenvalues of fij were 
negative, we would not be able to sample the operator as a Gaussian walk. To avoid this 
possibility, we implement the normal order process as follows: 

.3 



exp 



exp 



2V3 



T \M I exp 



e 
24 



(2 - (didjfij + di 



exp 



M | exp 



^ (~<W*,;) + y a ( 2 " v^) ( W« + a 



exp 



2y/S 



2^ 



T 



T 



(23) 



where A/" denotes the normal ordering of all derivative operators to the left. Since the 
left (and only the left) operator ex P(^g^) is already normal ordered with respect to the 
position-dependent operators in the middle term, the two normal ordered exponentials can 
be combined to remove the restriction of a positive definite fcj. Now, only the full covariance 
matrix C needs to be positive definite, which will always be the case for e sufficiently small. 
The final normal ordered exponential describes a non-uniform Gaussian random walk with 
mean /ij and covariance matrix Cjj : 



Hi 



a 



>■■:.) 



€ 

24 
e 

2^3 



1 



1 



(24) 
(25) 
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To sample this random distribution we need y/C, which we can approximate correctly to 
fourth order as 

1 ' 



i,3 



2y/S 



(26) 



Thus the entire factorization (|19D can be simulated by setting 

Wi = Xi + £ U 



l \2\ 

y l = w l (e/2)+i' l 

c 3 



1 



V3 r 



Zi = Vi 



e 
24 



. 2^3' 
(2-V3) Vl (y ) + ^ 



x[ = Zi(e/2) + 



where & to £f are four sets of independent Gaussian random numbers with zero mean and 
unit variance. 

As a severe test of the fourth order convergence of this algorithm, we use it to simulate 
the Brownian dynamics of 121 colloidal particles in two dimensions, with dimensionless 
surface density N/A = 0.5, interacting via a pairwise strongly repulsive Yukawa potential 

V(r) = —exp[-\(r-i)], (28) 
r 

with A = 8. This system has been described and simulated extensively via second order 
algorithms by Branka and Heyes []10f| . We will refer readers to this work for a detailed 
description of the system and their algorithms. In Fig. 1. we show the convergence of the 
potential energy at one parameter setting as a function of the time step-size used. (Compare 
this figure to that of Fig. 6 of Branka and Heyes ||10|| .) The linear and quadratic convergences 
are clearly evident. The two second order algorithms used are as described by (|ll|) and (|14|). 
These are referred to in Ref. as algorithms LGV2b and LGV2a respectively. 

When our fourth order Langevin algorithm is implemented by using the standard fourth 
order Runge-Kutta algorithm to solve the trajectory equation (Q) we obtained results as 
shown by open circles in Fig.l. The variance of the potential energy increases abruptly at 
around e = 0.0028 and the algorithm becomes unstable at larger e's. The problem can be 
traced to the instability of the Runge-Kutta algorithm itself in solving for the many-body 
dynamics. While the trajectory evolution exp(eD) should always decrease the potential 
energy, 

* = ar* = - |w| ' (29) 

this is no longer respected by the Runge-Kutta algorithm at larger time steps. The failure is 
due to the fact that Gaussian random walks can deposit particles so close together that the 
velocity field is changing too steeply for the Runge-Kutta algorithm to integrate accurately. 
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Each of these particles then gets placed chaotically somewhere in the periodic box, often 
again too near others, thus multiplying the number of particles that will be moved erratically 
in the next iteration. At time steps below but near e = 0.0028, the system can recover the 
regular behavior after several to hundreds of iterations, but only at the cost of increased 
variances and larger errors. Thus the inaccuracy in the trajectory determination causes the 
Langevin algorithm to fail prematurely. 

To improve on this situation, we monitor the difference between the results of the stan- 
dard fourth order Runge-Kutta and the embedded second order algorithm (|l^). We use the 
absolute value squared of this difference as a gauge of the fourth order method, even though 
it is strictly only an error estimate for the embedded second order algorithm. If the value 
of this difference is larger than some tolerance (0.01 in our case), we reject the result of 
the Runge-Kutta and recompute the trajectory more accurately by applying our trajectory 
algorithm twice at half the time step size. At small time steps, this incurs only a very small 
overhead. Even at a time step of 0.004 only 3% of the trajectories have to be re-evaluated. 
With this improvement, our fourth order Langevin algorithm gives results as shown by solid 
circles in Fig. 1. (We also applied similar monitoring processes to LGV2a and LGV2b by 
comparing the results of their first and second order Runge-Kutta algorithms.) The step- 
size dependence of the fourth order algorithm is remarkably flat, and yielded the converged 
results of the lower algorithms at step-sizes nearly 50 times as large. 

IV. SOLVING THE KRAMERS EQUATION 

While we are not aware of other multivariable 4th order Langevin algorithms, there 
are two fourth order algorithms in the literature for solving the Kramers equation in one 
dimension Ij7|,|nj|. Despite its more complicated appearance, the Kramers equation is actually 
simpler to solve than the Langevin equation. To illustrate the versatility of our operator 
approach, we will derive systematically a number of fourth order algorithms for solving this 
equation. Following Hershkovitz [71], we write the Kramers equations in the form 
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ft = F i (q)- 7 ft + Ci, (30) 

where the force is derivable from a potential, i^(q) = — <9jV^(q). A key simplification follows 
from the Hamilton form of the equation 

Qi =Pi 

pi = Fi(q) -jPi + Ci, (31) 
where Q is the zero-mean Gaussian random noise vector with variance 

(6(*)0(0> = |tM(*-0- (32) 

The advantage here is that the noise only affects the momentum, and classically, the mo- 
mentum commutes with the position-dependent force term. We will study the case of the 
bistable potential 

V(q) = q A - 2q\ (33) 

at parameter value 7 = 1 and (5 = 5. For each algorithm considered below, starting with 
q(0) = and p(0) = 0, we evolve the system to a finite time of t = 6. For comparison, we 
note that the total energy approaches the equilibrium limit of E = —0.8 at infinite time. 



Hershkovitz has formally derived a 4th order algorithm for solving ([Jl]) using Taylor 
expansion, but he has given an explicit implementation only for one dimension. In one 
dimension, each update of his algorithm requires one determination of the particle trajectory 
to 4th order, 4 Gaussian random variables, and one evaluation of the derivative of the force. 
The results of using his algorithm to evolve the system energy as a function of the time step 
size e is shown as solid squares in Fig.2. The standard 4th order Runge-Kutta algorithm, 
which requires four evaluations of the force, is used to solve for the particle's trajectory. 

To derive factorization algorithms in any dimension, we note that the probability density 
function evolves according to 

P(q,p,t) = LP(q,p,t), (34) 
10 



where 



7, 



L = jK + 7V P • p - p ■ V q - F(q) ■ V p = L x + L 2 + L 3 + L 4 . (35) 

To factorize the evolution operator exp(eL) for small e, we decompose L into exactly solvable 
parts T plus D and apply known fourth order factorization schemes |],[| . Drozdov and Brey 



TTH have recently initiated such a study of the Kramers equation. In this work, we have 
done an exhaustive search of all possible choices of solvable T and D such that [D, [T, D}] or 
[T, [D, T] is also solvable. We use the word "solvable" here loosely to denote either analytical 
result or trajectory determination. For example, the effect of exp[e(L2 + L3 + L4)] on the 
distribution function P(q, p, t) corresponds to evolving the particle trajectory forward in 
time with a linear friction. Since this can be computed using any trajectory integration 
algorithm, we consider L 2 + L 3 + L 4 to be solvable. While there are many solvable choices 
for T and D, such as the sum of any two Li, few resulting double commutators are simple. 
The possible choices for T and D are dramatically reduced if we insist that one of their 
double commutators is also structurally similar to the original T or D. There are then only 
three possibilities. 

The first possibility is to take 

T = L 1 + L 2 + L 3 , 

D = L 4 , (36) 



which is the choice originally made by Drozdov and Brey |TTJ. The Green's function corre- 
sponding to exp(eT) is known analytically (11 j, and can be sampled via 



Pi = 7e + [H, 

q , i = q i +Pi(l-e^ e )/'y + u i , (37) 

where corresponding to each pair of (pi, g»), (//j, Vj) is a pair of correlated Gaussian random 
numbers given by 
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Here, & and ^ are again two independent Gaussian random numbers with zero mean and 
unit variance. Note that at a given step size e, all the above functions involving e~ 7<E , etc., 
only need to be evaluated once at the beginning of the simulation. The operator exp(eD) 
can be exactly simulated by 

pi=p i + eF < (q). (39) 

As we will see, this choice is clever because there is no trajectory equation to solve. The 
double commutator required for a fourth order factorization is 

[D, [T, D]] = [L 4 , [L 3 , L 4 ]] = - V q |F| 2 • V p (40) 

which is just D but with a force V q |F| 2 . For each choice of T and D, there are three 
generic schemes [0] for factorizing the decomposed operator exp[e(T + D)} to fourth order 
with purely positive coefficients. For this choice of T and D, we found that schemes A and 
B of Ref. H give rather similar results, so we will only present results for schemes A and C. 
Scheme A and C are respectively, 

and 

e e{T+D) = ^T^eD^eT^eD^eT^eD^eT + Q ^ ^ 

where 

D = D + ^[D,[T,D]]. (43) 

The results of these two algorithms are shown as solid and open circles in Fig. 2. We will 
refer to these two as algorithms DB (Drozdov and Brey) and K4a respectively. Each al- 
gorithm evaluates the force three times and the derivative of the force once. Drozdov and 
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Brey's algorithm uses 4 Gaussian random numbers and K4a uses eight. For the extra effort, 
algorithm K4a has a much flatter convergence curve. Drozdov and Brey solved their one 
dimensional problem on a grid. We used Monte Carlo simulation, which can be generalized 
to any dimension. 

The second possibility is to take 

T = L 1 + L 2 , (44) 
D = L 3 + L 4 . (45) 

The operator exp(eT) now corresponds to an Ornstein-Uhlenbeck process in p, , 



p [= Pl e-^ + uUl-e-^), (46) 



P 

and exp(e.D) evolves the particle trajectory forward in time without friction, 



Pi =PiO) 



4i = *(e). (47) 

In this case, the simpler double commutator is 

[T,[D,T]] = [L 2 ,[D,L 2 ]] = - 1 2 D, (48) 

which does not require the derivative of the force. For this choice, we need to switch T «-> D 
in scheme A and slightly modify it as follows: 

e e{T+D) = e |eT e Ie[l-6 2 7 2 /72]Z) e § e r e I e [l- e V/72]O e ieT + ^ 

The effect of the double commutator simply reduces the time of the trajectory evolution. 
This algorithm, which will be referred to as K4b, requires two trajectory determinations but 
no derivative of the force and only three Gaussian random numbers. The trajectory can be 
computed using the standard 4th order Runge-Kutta algorithm with four force evaluations, 
or the 4th order Forest-Ruth symplectic algorithm |T3[ with three force evaluations. The 



results from these two cases are plotted as solid and open diamonds respectively in Fig.2. 
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For this choice of D, we did not bother with factorization schemes B or C, since either would 
have required more than two trajectory determinations. 
The third possibility is to take 

T = L 1 , (50) 
D = L 2 + L 3 + L 4 , (51) 

where now exp(eT) is just a Gaussian process in pi, 

p'i=Pi + C*V^, (52) 

and exp(e.D) evolves the particle trajectory forward in time with friction. For this case, we 
have the simplest result, 

[T,[D,T\]=0, (53) 
and a simplified fourth order factorization 

e e(T+D) = ^T^eD^eT^eD^eT + (54) 

We shall refer to this as algorithm K4c. This algorithm is similar to K4b, with no force 
derivative necessary. If we solve the trajectory equation by the 4th order Runge-Kutta 
algorithm, we obtain results as shown by solid triangles in Fig.2. Note that in contrast 
to previous algorithms, this algorithm does not converge monotonically. It overshoots and 
converges from the top. 

In the course of our calculations, we find that for each algorithm, a more accurately 
determined particle trajectory will yield a flatter convergence curve. If we now further 
decompose D = D\ + D 2 in algorithm K4c, with 

A = L 2 , 

D 2 = L 3 + L 4 , (55) 
the double commutator [D ly [D 2 , Di]] = — r y 2 D 2 is just a restatement of (H3). We can again 



factorize, 

14 



The friction evolution e eDl rescales the momentum, 

Pi = Pi^% (57) 

and e eD2 again evolves the trajectory forward for time e. This way of solving the trajectory 
with friction doubles the number of trajectory calculations, but also further flattens the 
convergence curve. To minimize the number of force evaluations, we use the Forest-Ruth 
symplectic algorithm to calculate the trajectory. The results are shown as open triangles in 
Fig.2. 

Of the algorithms studied, Drozdov and Brey's algorithm makes maximum use of ana- 
lytical knowledge and is very efficient. The improvement we suggested, algorithm K4a, with 
twice the number of Gaussian random numbers, seemed to double the range of the conver- 
gence. Our new algorithms K4b and K4c, while requiring two trajectory determinations, 
have no need of evaluating the force derivative. All these algorithms serve to illustrate the 
power of the factorization method. While the diligence of Hershkovitz is rewarded with just 
a single fourth order algorithm, we can survey the form of the evolution operator and derive 
many fourth order algorithms. 

V. SUMMARY AND CONCLUSIONS 

In this work, we have shown how the method of operator factorization can be applied to 
the Langevin equation to derive a practical fourth order algorithm. This method of factor- 
izing an evolution operator of the form e e ( A+B ) leads to unitary algorithms for solving the 
Schrodinger equation in quantum mechanics, symplectic algorithms for solving Hamilton's 
equations in classical mechanics, and norm preserving algorithms for solving the Langevin 
equation in stochastical mechanics. A key step in deriving a fourth order Langevin algorithm 
is our treatment of the double commutator term through successive use of normal ordering. 
The resulting algorithm (|27|) is computationally demanding, but one is rewarded by a very 
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flat convergence curve, virtually eliminating the step-size dependent error. Future use of 
this algorithm in other applications may lead to further simplifications and enhancements 
of its utility. 

We also derived a number of 4th order algorithms for solving the Kramers equation. 
The freedom in decomposing the kernel operator and choosing a particular factorization 
scheme illustrates the power of this approach. It is difficult to see these global structures 
from just doing Taylor expansions. One advantage of our simulation approach is that we 
are not restricted to solving the Kramers equation in one dimension. We can solve it in any 
dimension. Our use of the Kramers equation is also only illustrative, one can apply this 
method of operator factorization to other stochastic equations of one's own interest. 

It is observed in solving both equations that the step-size error is reduced by solving 
the trajectory more exactly. Different fourth order algorithms for solving the trajectory 
equation can yield different convergence curves. One should therefore explore the effect of 
using fourth order algorithms other than Runge-Kutta in implementing any of the above 
stochastic algorithms. 

ACKNOWLEDGMENTS 

This research was funded, in part, by the U. S. National Science Foundation grants 
PHY-9512428, PHY-9870054 and DMR-9509743. 



16 



REFERENCES 

[1] H. Risken, The Fokker-Planck Equation, Methods of Solution and Applications, 2nd 
edition, (Springer, New York, 1989). 

[2] E. Helfand. Bell Syst. Tech. J., 58, 2289 (1979). 

[3] I. Drummond et. al, Nucl. Phys. B220, 119 (1983). 

[4] A. Ukawa and M. Fukugita, Phys. Rev. Lett. 55, 1854 (1985). 

[5] S. A. Chin, Nucl. Phys. B (Proc. Suppl.) 9, 498 (1989). 

[6] S. A. Chin, Phys. Rev. A42, 6991 (1990). 

[7] E. Hershkovitz, J. Chem. Phys. 108, 9253 (1998). 

[8] M. Suzuki, Computer Simulation Studies in Condensed Matter Physics VIII, eds, D. 
Landau, K. Mon and H. Shuttler (Springier, Berlin, 1996). 

[9] S. A. Chin, Phys. Lett. A226, 344 (1997) . 

[10] A. Branka and D. Heyes Phys. Rev. E60, 2381 (1999) . 

[11] A. N. Drozdov and J. J. Brey, Phys. Rev. E57, 1284 (1998) . 

[12] M. Suzuki, J. Math. Phys. 32, 400 (1991). 

[13] E. Forest and R. D. Ruth, Physica D 43(1990) 105. 



17 



FIGURES 
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FIG. 1. The convergence of Langevin algorithms for simulating the Brownian dynamics of 121 
interacting colloidal particles in two dimensions. The equilibrium potential energy per particle 
is plotted as a function of the time step size e used. Open diamonds are results using the first 
order Langevin algorithm. Solid triangles and solid squares denote results of the two second 
order algorithms LGV2a and LGV2b respectively, as described in the text. Open circles give 
results of our fourth order Langevin algorithm using the standard 4th order Runge-Kutta algorithm 
for determining the particle trajectory. The solid circles give results with improved trajectory 
determination as discussed in the text. 
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FIG. 2. The convergence of various fourth order algorithms for solving the Kramers equation 
in one dimension. The energy calculated is at a finite time of t = 6 with system parameters (3 = 5 
and 7 = 1. Solid squares: Hershkovitz's algorithm. Solid and open circles: Drozdov and Brey's 
algorithm and K4a. Solid and open diamonds: two variants of algorithm K4b. Solid and open 
triangles: two variants of algorithm K4c. See text for algorithm descriptions. The fitted lines all 
have leading term e 4 or higher. Error bars are comparable or smaller than the size of plotting 
symbols. 
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